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We report on automating the Catani-Seymour dipole subtraction which is a general procedure to treat infrared divergences 
in real emission processes at next-to-leading order in QCD. The automatization rests on three essential steps: the creation 
of the dipole terms, the calculation of the color linked squared Born matrix elements, and the evaluation of different helicity 
amplitudes. The routines have been tested for a number of complex processes, such as the real emission process gg — > ttggg. 



1. INTRODUCTION 

QCD as the gauge theory of the strong interaction 
allows to predict cross sections for hard scattering re- 
actions which, at a hadron collider, typically involve 
high multiplicities of colored partons. In the perturba- 
tive approach calculations based on exact QCD matrix 
elements at leading order (LO) provide first estimates 
for cross sections and differential distributions. How- 
ever, in the complicated environment of a hadron col- 
lider precision calculations to next-to-leading-order 
(NLO) in QCD are often needed in order to reli- 
ably predict (and separate) the Standard Model back- 
ground from possible new physics signals. In the era 
of LHC this has triggered a lot of activity concerning 
NLO QCD corrections to multi-particle reactions, see 
e.g. Refs. mi. 

The salient feature of NLO corrections is the 
presence of virtual and real emission contributions. 
Virtual-loop corrections exhibit both ultraviolet (UV) 
and soft and collinear divergences which we call IR 
divergencies in the following. The real corrections 
contain only IR divergencies from soft and collinear 
emissions. Upon summation of the two parts all 
IR divergences cancel (for so-called IR safe observ- 
ables) 03 1415 1 . Since virtual and real corrections have 
different phase space integrals, these cancellations are 
not always trivial. In the Catani-Seymour dipole for- 
malism Il6l7l8l . the IR divergences of virtual and real 
corrections are treated separately by subtracting suit- 
able dipole terms so that each of the contributions be- 
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comes individually finite. The dipole terms are con- 
structed systematically relying on the universal nature 
of soft and collinear limits in QCD. Thus the method 
allows for a general treatment of IR divergences to 
NLO in QCD. 

Current applications in phenomenology consider 
processes with six or more parton legs ri'21 which 
require about one hundred dipole terms. These cal- 
culations are rather tedious and since the algorithm 
underlying the Catani-Seymour dipole subtraction is 
a combinatorial one automatization is favored. The 
construction of the complete subtraction terms re- 
lies to a large extent on squaring color correlated 
Born amplitudes and dressing them with the corre- 
sponding dipoles. To that end, we can use existing 
software for the automatic evaluation of Born ampli- 
tudes 1 9 10 11 12 13 14 15 161 by means of a suitable 
interface. Thus, the completely automatic generation 
of all subtraction terms in the dipole formalism be- 
comes feasible. 

In related work Ref. IfTTl recently reported on de- 
tails of an automatization, although the code is un- 
published. Ref. [18J made code publicly available. 
However, all process dependent information, i.e. the 
(color correlated) Born squared matrix elements, still 
have to be provided by the user. Automating this step 
is precisely what we are aiming at in the present arti- 
cle. 

2. ALGORITHM 

In this section we briefly review the algorithm of 
the dipole subtraction with particular emphasis on the 



features of the real emission contributions for a given 
process. 

1. Choose all possible emitter pairs from the ex- 
ternal legs. In the dipole subtraction, the root of 
the splitting of the quarks and gluons is called emit- 
ter For convenience we call the two fields into 
which an emitter splits, emitter pair We use in- 
dices, ij, and k, for fields in a final state and, re- 
spectively, indices, a and b, for an initial state field. 
The quark (anti-quark) and gluon are denoted by / 
(/) and g. In case both partons of an emitter pair are 
in the final state, possible combinations are (/, /) = 
(!)(/,§), (2)fe,g),(3)(/,/). In case of one paiton 
in the initial and the other in the final state, we have 
(fl,/) = {4)if,g),{5)ig,g)A6){fJ),{l){gJ). There 
are also the other combinations where the quarks are 
replaced by the anti-quarks in the cases, (1), (4), (6), 
and (7). 

2. Choose all possible spectators for each emitter 
pair The spectator is one external field which is dif- 
ferent from both fields of the emitter pair For a spec- 
tator in the final (initial) state denoted by k (b), this 
condition means k ^ i,j (b ^ a). It emerges from a 
special feature of the subtraction formalism namely 
that the color factors of the square terms |M,p are 
expressed through the ones of the interference terms 
M,My (/ 7^ j) due to color conservation. 

3. Construct the dipole terms from the chosen com- 
binations of emitter and spectator. The previous steps 
provide all such combinations as pairings (emitter, 
spectator)= {ij,k),{ij,b),{ai,k), and {ai,b). Each 
case corresponds to one dipole term, D,y ,t,D" ,D^', 
and D"'*, respectively, and explicit expressions are 
given in [6*81. For example, the dipole term D,yjt in 
the massless case reads, 

^ -1 ... ,_,Ta.-T, 



2pi-Pj 
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(1) 



The quantity, {ij^k JTa- -T,,! ij^k), is called color 
linked Born squared matrix element (CLBS). It is 
given by the Born amplitude squared with two ad- 
ditional color operator insertions at the emitter and 
spectator legs. The color operator T denotes either 
a fundamental f-' or an adjoint /"*' generator, de- 
pending on the parton type (i.e. quark, anti-quark or 
gluon). The quantity ^ij.k is the so-called dipole split- 
ting function. In case of the emitter being a quark, for 
example, [ij^k) = {fg,k) the splitting function is di- 



agonal in the spin space of the quark: 

(^|V/,i|.^') = 
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where z, and ytji^ are functions of the external mo- 
menta as. 



Zi = 



Pi ■ Pk 



yij.k 



Pj ■ Pk + Pi ■ Pk 
Pi ■ Pi 



Pi-Pj+Pj-Pk+Pk-Pi 



(3) 
(4) 



In case of the emitter being a gluon, for example, 
{ij,k) — {ff,k), it does exhibit a correlation with the 
gluon helicity according to 



(A'|V//,A.|v) = 87trfia.v 



r,/JV 
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(5) 



which has to be treated accordingly (see for example 
Ref. 1 191). 

In summary, these three steps generate all dipole 
terms to subtract all IR divergences of a real emission 
process at NLO. The construction of the integrated 
dipole terms is less involved and proceeds in complete 
analogy 1,6,8] , thus we skip a detailed discussion here. 

3. CODE STRUCTURE 

Here we briefly sketch our implementation of the 
dipole subtraction formalism. There is a large free- 
dom how such an implementation can be done. In 
Ref. [201 for example the implementation was done 
in form of two independent C/C-H- libraries provid- 
ing all the necessary functions to evaluate the dipole 
terms. The slight disadvantage of this approach is that 
the produced code is non-local and that there is some 
redundancy in the calculation. In the present work 
we follow a different approach. The main idea here 
is to have code generator which will produce an op- 
timized flat code which can be further optimized by 
the compiler To do so we haven chosen to interface a 
Mathematica program with MadGraph li9ilOJ . 

3.1. Mathematica code 

We implement the creation of the dipole terms in 
Mathematica. With a given (real emission) parton 



scattering process as an input, the Mathematica code 
automatically writes down all dipole terms needed at 
NLO. It provides all expressions explicitly except for 
the CLBS. The code creates the dipole terms in an or- 
der according to the kind of the emitter pairs, i.e. the 
seven combinations of {i,j) or {a,i) listed in Sec.|2] 

The first group of dipole terms (dipole 1 ) are the 
ones with the emitter pairs, (1),(2),(4), and (5). These 
emitters reduce the NLO real emission process to a 
Born amplitude which is the LO contribution to a pro- 
cess with one less gluon in the final state. The second 
kind of dipole terms (dipole 2) has the emitter pair 
(3), while the third and fourth (dipole 3 and 4) have 
(6) and (7) as emitter pairs. 

In order to demonstrate the code, let us discuss the 
example g{a)g{b) -^ u{\)u{2)g{'i). The code starts 
with the creation of the first dipole Dn 2 which be- 
longs to the group dipole 1 and the output is written 
in the form. 



D,;^/gi(l32) 



-1 



-Vow(132) 



Bl(132) 



T2 



(6) 



2p; • Pi 

where Vijkfgk is the dipole splitting function and B 1 
denotes the CLBS. The indices, ijkfgk{l32), of D 
and V mean that (emitter, spectator)=(/y,A;)=(quark 
gluon, something)=(13,2). About the CLBS the code 
stores only the necessary information for the direct 
calculation. It writes each CLBS as B'i' correspond- 
ing to Dipole 'i', where i can be 1,2,3, or 4. For in- 
stance, the output for Bl(132) in Eq. Q is returned 
in the form, 

Bl(132)=Bl[{{g,pa},{g,pb}}--> 

{{u,pijtil[l,3]}, {ubar,pktil[2]}}] , (7) 

where gg — > uu is the reduced Born process. The 
function pijtil[l,3] is the reduced momenta for the 
emitter and pktil[2] for the spectator In general, for a 
given NLO real emission process with n parton legs, 
each dipole term has a reduced Born squared matrix 
element with (n — 1 ) parton legs. The reduced (n — 1 ) 
external momenta are functions of the original n ex- 
ternal momenta. For example, the reduced momen- 
tum for an emitter in the dipole term D,j ,t reads 






^r 



-pI 



(8) 



The Mathematica code provides explicit expressions 
for the reduced momenta of each dipole in the output 
according to Eq. (|7]). 



At the end of the run, the total number of generated 
dipoles is shown. The final output for the complete 
subtraction term is given as a C- or Fortran routine for 
numerical evaluation. We have tested the generation 
of dipoles with our Mathematica code for complex 
processes, obtaining, e.g. twenty seven dipole terms 
for the process gg -^ uug, eighty dipoles for gg -^ 
uuddg, and one hundred dipoles for gg -^ uuggg. 

3.2. Interface to MadGraph 

As we see, the remaining ingredient at this stage is 
the CLBS appearing in all dipole terms. In case the 
emitter is a gluon, they also include the different he- 
licity components of the CLBS for the emitter. In or- 
der to obtain these quantities in an automatic way, it is 
advantageous to use a publicly available software for 
automated LO calculations. We choose MadGraph 
for this purpose and interface our Mathematica pro- 
gram with the stand-alone version 1.9. lOJ . 

Let us briefly explain our interface to MadGraph to 
obtain the CLBS. In MadGraph, the color factors are 
separated from each diagram. In the evaluation every- 
thing is expressed in terms of generators of the funda- 
mental representation. A typical example is that the 
factor /"'"^ of the gluon three point vertex is rewritten 
in terms of the fundamental generator f ■' due to the 
identity, 

jrabc ^ _2/( Tv[t"t''t'] - Tr\fth''\ ) . (9) 

The color factors of each diagram are sorted in an 
unique order and they are expressed in a sum of some 
terms. When a specific term of a diagram is identical 
to one of the other diagrams, it is combined as 



M = ^Q,J„, 



(10) 



where Ca denotes the independent color factors. Each 
Ca has fundamental and adjoint color indices corre- 
sponding to the external quarks (anti-quarks) and the 
gluons, respectively. !„ is the joint amplitude, e.g. 
Ji = +Ai — A3 + • • • where A,- is the partial ampli- 
tude of /-th diagram (with the color factor stripped 
off). The invariant matrix element squared is finally 
expressed in the form. 



IMI- 



J)^CFJ, 



where the color matrix CF is defined as 
(CF)„, = ^ CQ. 



(11) 



(12) 



For the CLBS we need to evaluate Eq. ( fTTI ) with 
an insertion of two additional color operators to the 
emitter and spectator legs. This is precisely what our 
interface to MadGraph does. 

Let us return to the example of B 1(123) from 
the previous subsection. The reduced Born pro- 
cess g{a)g{b) — * m(1)m(2) has three diagrams and 
the color factors are combined into two indepen- 
dent ones, (CuCi) = ((f"f'')i2,(f''f")i2)- The com- 
ponents of the color matrix are written in the 
traces, (CF)ii = (CF)22 = Tr[rVf"f''] and (CF)i2 = 
(CF)2i = TT[t''t"t''t"]. Then the color matrix is calcu- 
lated as 

16/3 -2/3 
-2/3 16/3 

For the CLBS Bl(123) we need the fundamental op- 
erator insertions into the legs 1 and 2. The compo- 
nents of the color matrix CF are modified to (CF') 1 1 = 
Tr[t''t"ft"fy''] and (CF')i2 ^ Tr[t''t"ffy"t% Then 
the modified color matrix is calculated as 



CF 



(13) 



' 10/9 1/9 



(14) 



The subroutines of MadGraph for the color factor cal- 
culations are well structured and the original routines 
to add the color factors tf- and /"*' can be applied to 
the additional color insertions for CLBS. We have re- 
alized the two color insertions in an automatic way 
and checked that MadGraph with our interface works 
for rather involved processes. One of the most com- 
plex checks consists of the two color insertions into 
the process g{a)g{b) -^ u{\)u{2)g{3)g{A). In Mad- 
Graph the normal color matrix for the process is a 24 
by 24 matrix. Here we show only the first 15 compo- 
nents in the first row as 

CF=— (512,8,-64,80,8,-10, (15) 

-1,-64, -64,8, -1,-10, -1,62, -10, •••). 

Next, we perform two adjoint operator insertions into 
the legs 3 and 4, and the extended routines calculate 
the modified color matrix as 



CF' = -(8,0,8, 16,0, -2, 



(16) 



0,8, -1,-1, 1,2, -8, -7, 1,---). 

We have checked that the result in Eq. ( fT6] l and the 
sum of all components agree with results of our inde- 
pendent private code. 



Next we briefly comment on the different helicity 
components of the CLBS for the gluon emitter In 
MadGraph the gluon polarization vector is calculated 
by the subroutine ' VXXXXX' of the HELAS library 
where the polarization vector is taken in the circular 
polarization representation 112 11221 . Then, it is favor- 
able to calculate the dipole terms with the helicity cor- 
relation in the circular polarization basis. For exam- 
ple, in the previous process g{a)g{b) -^ u{\)u{2)g{'i), 
we take one dipole D"^ and calculate it in the circu- 
lar polarization basis as V"^ (A,, A-')B 1 (g (X, A-')g — * mm) 
where V'i'^(A-,A,') is constructed from the dipole spUt- 
ting function in the basis of the Lorentz indices as in 
Eq. (|5]l by multiplying the circular polarization vec- 
tors. The arguments (X,A-') are the ones for the dif- 
ferent helicities of the original and the complex con- 
jugated amplitude for the emitter gluon. We have 
completed the interface to obtain the CLBS with the 
different helicities as Bi(A-,A,') and checked it in pro- 
cesses as the above Bl(g(A,,A-')g -^ uii). 

3.3. Complete structure 

Let us finally wrap up by displaying the complete 
structure of our code shown in the flowchart in Fig. [T] 
An user specifies a n-parton real emission process 
as an input to the Mathematica code as well as a 
set of the external momenta {pn}- The Mathemat- 
ica program then generates all dipole terms in the ap- 
propriate order (see Sec. 13. Il l along with the CLBS 
B 1 , ■ • • , B4. The created dipole terms are written as 
C- or Fortran routines to a file 'dipole. c' or 'dipole.f. 
In this file the explicit expressions for the dipole split- 
ting function V are contained. Together with each 
dipole, the information on its reduced kinematics as 
a function of the external momenta is stored (see 
Sec.imi. 

The necessary information to calculate each 
CLBS is transferred to MadGraph through the file 
'proc_card_dip.data'. This data file is an imitation 
of the input file 'proc_card.data' of the original Mad- 
Graph. Our interface reads 'proc_card_dip.data' and 
gets MadGraph to write Fortran routines for the eval- 
uation of each CLBS in a file 'matrix_dip.f ' . As we 
see in Sec. 13.21 the color matrix CF is modified to CF' 
due to the additional color insertions. 

Finally, the C- and/or Fortran codes in the two files 
can be used for the numerical evaluation of all dipole 
terms as functions of the external momenta {p,,}. The 



Input[gg ^uiig,{p5}] 



Mathematica code 



pioc_caid_dip.data 



dipole.f 



D(1) = V(132)B1{132) 

Bl(132) = Bl[gg ^ 11(13) u(2)] 
{P4> = fi({P5}) 



D(13)-V(123,A,X')B2(123.X,X') 
B2(123) = B2[gg^g(]2,;t,X')g(3)]- 



D(27) = V(b2a)B4(b2a) 
B4Cb2a) = B4[ g(a) vfb2) ^ u g ] 



Interface to Madgraph 



matrix_dip.f 



Bl(132)= Jbi CFi' Jbi 



>B2(123) =7b2(X,')CFi'3Tb2W 



B4(b2a)=TBlcF27lB4 



Figure 1. 

The complete structure of the code is shown. The Mathematica generates all subtraction terms in the dipole 
formalism and returns output in C or Fortran for numerical evaluation. The calculation of the CLBS is delegated 
to MadGraph via an interface. 



sum of the dipole subtractions from the invariant ma- 
trix element squared of the NLO real emission pro- 
cess reads. 



|M(2 ^ (n - 2)partons)p - ^D(0 . 



(17) 



The invariant matrix element squared can be calcu- 
lated with the original MadGraph version and / runs 
over all dipole terms. Eq. ( [TtT i is finite upon integra- 
tion over the phase space of the unresolved parton. 
Thus, it can finally be integrated over the phase space 
by using standard Monte Carlo techniques to obtain 
its contribution to an IR safe cross section. 



4. OUTLOOK 

We have reported on ongoing work to automate the 
Catani-Seymour dipole formalism in order to calcu- 
late the subtracted invariant matrix element squared 
Eq. ( [TtI i in an automatic way. 

The automatization essentially requires three ingre- 
dients: the automatic generation of all dipole terms, 
the calculation of the CLBS, and the evaluation of 
different helicity amplitudes. The implementation of 
each of these tasks either in our Mathematica program 
or in an interface to MadGraph has been completed 
and the respective routines have undergone sufficient 
checks. We are now finalizing the user interface and 
the output format to achieve full automatization. At 



the same time we are checking our code for vari- 
ous massless real emission processes, like gg -^ uug, 
gg — > uugg, and gg -^ uuggg, as well as for massive 
processes, hke gg -> ttg, gg -^ ttgg, and gg -^ ttggg 
to obtain finite results for |Mp — ^D in all soft and 
collinear limits (and IR safe contributions to cross 
sections). Once the reliability of our software has 
been fully established and the code has been opti- 
mized for speed, it will be made publicly available. 
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